home *** CD-ROM | disk | FTP | other *** search
- /****** Moon.rexx *************************************************************
- *
- * $VER: Moon 2.0 (30.4.95)
- *
- *******************************************************************************
- *
- * Originially by Greg Cunningham 1990
- * Rewritten by kas@mink.ping.dk (Klaus Alexander Seistrup) 1995
- *
- * Converted from Unix TODAY.C
- *
- * GMT settings (America):
- *
- * Eastern 5
- * Central 6
- * Mountain 7
- * Pacific 8
- *
- * Denmark 1
- *
- ******************************************************************************/
-
- TimeZone = 1 /* Central European Time */
- DaylightSavings = 1 /* Set to 1 if daylight savings time */
-
- rmLib = 'rexxmathlib.library'
-
- IF ~SHOW(L,rmLib) THEN CALL ADDLIB(rmLib,0,-30,0)
-
- SAY Lunar(DATE('S'),TIME(),TimeZone-DaylightSavings)
-
- EXIT 0
-
- /*
- ** Phase of the Moon. Calculates the current phase of the moon.
- ** Based on routines from `Practical Astronomy with Your Calculator',
- ** by Duffett-Smith.
- **
- ** SYNOPSIS:
- **
- ** string = Lunar(year,time,gmt)
- ** date -- Current sorted date
- ** time -- Current time of the day
- ** gmt -- Hours from GMT time zone (subtract one if daylight savings)
- **
- ** say Lunar(DATE('S'),TIME(),0) Current phase of the moon
- ** say Lunar(DATE('S')+1,TIME(),0) Same time tommorow
- ** say Lunar('19900225','12:00:00',0) Feb 25, 1990 at noon
- */
- Lunar: PROCEDURE
- ARG cdate,ctime,GMToffset
- /*
- ** Days between 01-Jan-1985 and now
- */
- numDays = DATE('B',cdate,'S') - DATE('B','19850101','S') + 1
- /*
- ** Add fraction of day
- */
- PARSE VAR ctime hour':'min':'sec
- hour = hour + GMToffset
- IF (hour > 23) THEN
- DO
- numDays = numDays + 1
- hour = hour - 24
- END
- /*
- ** Add on percentage of day
- */
- numDays = numDays + ((hour + (min / 60) + (sec / 3600)) / 24)
- phase = PotM(numDays)
- cp = 'The moon is'
- IF (phase < 3.125) THEN cp = cp 'new'
- ELSE IF (phase >= 96.875) THEN cp = cp 'full'
- ELSE IF (phase >= 46.875) & (phase < 53.125) THEN
- DO
- phase2 = PotM(numDays + 1)
- IF (phase2 > phase) THEN cp = cp 'at the first quarter'
- ELSE cp = cp 'at the last quarter'
- END
- ELSE IF (phase >= 53.125) THEN
- DO
- phase2 = PotM(numDays + 1)
- IF (phase2 > phase) THEN cp = cp 'waxing'
- ELSE cp = cp 'waning'
- cp = cp 'gibbous (' || TRUNC(phase) || '% full)'
- END
- ELSE
- DO
- phase2 = PotM(numDays + 1)
- IF (phase2 > phase) THEN cp = cp 'waxing'
- ELSE cp = cp 'waning'
- cp = cp 'crescent (' || TRUNC(phase) || '% full)'
- END
- RETURN cp
-
- /*
- ** Calc phase-of-the-moon; returns percentage of full
- */
- PotM: PROCEDURE
- ARG numDays
- EPSILONg = 279.611371 /* Solar ecliptic long at EPOCH */
- RHOg = 282.680403 /* Solar ecliptic long of perigee at EPOCH */
- e = 0.01671542 /* Solar orbit eccentricity */
- Lzero = 18.251907 /* Lunar mean long at EPOCH */
- Pzero = 192.917585 /* Lunar mean long of perigee at EPOCH */
- Nzero = 55.204723 /* Lunar mean long of node at EPOCH */
- PI = 3.1415926535897932384626433
- N = 360 * numDays / 365.24219878 /* Days*360/Real_Days_of_a_Year */
- N = Pos360(N) /* Round up to positive 360 degrees */
- Msol = N + EPSILONg - RHOg
- Msol = Pos360(Msol)
- Ec = 360 / PI * e * SIN(Deg2Rad(Msol))
- LambdaSol = N + Ec + EPSILONg
- LambdaSol = Pos360(LambdaSol)
- l = 13.1763966 * numDays + Lzero
- l = Pos360(l)
- Mm = l - (0.1114041 * numDays) - Pzero
- Mm = Pos360(Mm)
- Nm = Nzero - (0.0529539 * numDays)
- Nm = Pos360(Nm)
- Ev = 1.2739 * SIN(Deg2Rad(2 * (l - LambdaSol) - Mm))
- Ac = 0.1858 * SIN(Deg2Rad(Msol))
- A3 = 0.37 * SIN(Deg2Rad(Msol))
- Mmprime = Mm + Ev - Ac - A3
- Ec = 6.2886 * SIN(Deg2Rad(Mmprime))
- A4 = 0.214 * SIN(Deg2Rad(2 * Mmprime))
- lprime = l + Ev + Ec - Ac + A4
- V = 0.6583 * SIN(Deg2Rad(2 * (lprime - LambdaSol)))
- ldprime = lprime + V
- D = ldprime - LambdaSol
- perc = 50 * (1 - COS(Deg2Rad(D)))
- RETURN TRUNC(perc+0.005,2) /* Return percent rounded to two decimals */
-
- /*
- ** Get a positive 360 degree mod
- */
- Pos360: ARG deg
- deg = deg // 360
- IF (deg < 0) THEN deg = deg + 360
- RETURN deg
-
- /*
- ** Convert degrees to radians
- */
- Deg2Rad: PROCEDURE
- ARG deg
- PI = 3.1415926535897932384626433
- RETURN (deg * PI) / 180
-
- /*
- ** EOF
- */
-